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I.   INTRODUCTION  AND  SUMMARY 

A.   INTRODUCTION 

The  failure  rate  function,  or  hazard  function  (hazard 
for  short)  may  be  described  as  the  conditional  probability 
of  an  equipment's  failing  at  operating  age   t,  having  sur- 
vived to  that  age.   The  reliabilities  of  a  variety  of  elec- 
tronic and  mechanical  items  are  conveniently  and  naturally 
described  in  terms  of  the  appropriate  hazard  function,  and 
so  is  the  longevity  of  human  beings.   The  term  force  of 
mortality  replaces  hazard  in  the  latter  context. 

This  paper  is  devoted  to  a  study  of  several  simple 
analytical  representations  for  hazard  functions.   These 
representations  are  in  turn  based  upon  representations  of 
random  variables  having  certain  required  properties,  in 
terms  of  others  having  familiar  distributions — in  particular 
the  exponential.   Similar  ideas  are  due  to  Tukey  [1],  and 
recently  have  been  examined  by  Parzen  [2].   The  hazard 
representations  proposed  are  quite  expeditiously  used  in 
simulation  studies,  e.g.  of  system  reliability  or  avail- 
ability in  terms  of  component  lifetimes.   They  may  also  be 
used  in  data  analysis  studies,  in  order  to  parsimoniously 
describe  data  sets  in  terms  of  perturbations  of  convenient 
and  familiar  standard  distributions.   Their  use  in  data 
analysis  and  simulation  is  also  described  in  Gaver,  Laven- 
berg,  and  Price  [3],  and  in  Gaver  and  Chu  [4]. 


B.   SYSTEM  FAILURE  PATTERNS 

It  is  plausible  to  think  that  the  time  series  of  fail- 
ures in  a  system  may  involve  these  stages. 

Early  failures.   There  may  be  a  relatively  large  number 
of  failures  soon  after  a  system  is  introduced  because  of 
design  defects,  production  errors,  or  errors  stemming  from 
maintenance  personnel  inexperience.   This  situation  is 
characterized  by  a  hazard  function  that  is  initially  large, 
but  that  decreases  with  time.   "Infant  mortality"  is  in 
evidence. 

Random  Failures.   Following  the  early  failure  period 
there  may  be  a  period  during  which  failures  occur  at  an 
essentially  constant  rate  for  a  rather  prolonged  time. 
During  this  period  the  hazard  function  is  nearly  constant, 
so  the  times  between  failures  are  close  to  being  exponen- 
tially distributed.   The  effect  of  age  or  wearout  is  not 
yet  apparent. 

Wearout  Failures.   Eventually  following  the  period 
during  which  a  constant  hazard  is  evident  there  is  likely 
to  be  a  period  of  ever-increasing  failure  rate  caused  by 
wearout  of  system  components. 

A  graphical  representation  of  a  hazard  function  that 
exhibits  the  behavior  described  is  given  below.   Note  that 
it  has  the  legendary  "bath  tub"  shape. 


> 

U 

3 

u 

(0 

-Q 
Eh 

x: 

-p 
(d 

CQ 


•H 

fa 


43 


Some  comments  on  the  above  follow: 
The  term  "failure"  may  refer  to  an  event  that  is 
analogous  to  human  death,  after  which  the  entire  system 
is  replaced.   On  the  other  hand  repair  or  component  replace- 
ment may  occur  after  failure:   the  system  is  only  repaired, 
not  entirely  replaced.   In  the  former  case,  a  hazard  func- 
tion of  the  kind  depicted  in  Figure  1  applies  to  each  system 
event  ("death");  when  the  system  is  installed  (or  is  born), 
that  hazard  operates  starting  from  scratch  at   t  =  0   until 
system  failure  (human  death,  for  instance) ,  after  which  a 
similar  hazard  goes  into  effect,  starting  once  again  from 
zero.   In  the  latter  case,  in  which  repair  of  a  component 
occurs,  a  hazard  function  like  that  of  Figure  1  applies  at 
t  =  0,  but  after  the  first  event  ("failure")  at   t1   a  re- 
pair action  is  accomplished.   The  same  hazard  operates  for 
t  _>  t,   until  the  next  event  at   t2  >  t,  ,  and  so  on.   Inter- 
mediate situations  may  be  envisioned,  in  which  after  event 
n   at   t   the  hazard  governing  system  failure   n+1   starts 

at   t   -t,0<t   <t. 
n    n       n    n 

Although  there  is  reason  to  assume  that  hazards  some- 
what like  that  of  Figure  1  occur  in  general  for  systems , 
the  possibility  exists  that  the  system  hazard  is  "bumpy" 
because  wearout  failures  of  components  or  subsystems  may 
well  occur  at  intermediate  times. 

If  the  theory  is  applied  to  systems  with  little  or  no 
wearout  propensity,  as  should  be  the  case  when  dealing  with 
computer  software  modules,  then  the  hazard  function  may 
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well  exhibit  the  initial  falloff  of  Figure  1  but  not  the 
rise  at  later  times.   In  fact,  a  constant  decline  as  bugs 
are  found  and  removed  could  be  (optimistically)  anticipated 
for  software.   The  right-hand  side  of  the  bath  tub  vanishes, 
and  the  picture  is  that  of  a  ski  slope. 
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II.   ANALYTICAL  HAZARD  REPRESENTATIONS 

A.   MODELS  FOR  THE  HAZARD  FUNCTION 

In  this  section  mathematical  models  are  presented  for 
the  failure  rate  or  hazard  function.   Recall  that  the  hazard 
may  be  defined  as  follows. 

Definition.   Suppose  that  the  time  to  failure,  X,  is  a 
random  variable  with  distribution  function  F(x),  where 
F(0)  =  0;  the  latter  possesses  the  density  function   f(x), 
f(x)  =  dF/dx,  such  that  for  any  positive   x, 


x 

F(x)  =  /   f(y)  dy  .  (2.1) 

0 


Then  the  hazard  function ,  or  failure  rate  at  age   x,  is 
given  by 

h(x)  -  1  -(F(x>  •  (2-2) 

The  interpretation  of   h(x)dx   is  that  it  is  the  con- 
ditional probability  of  failure  in  the  interval   (x,  x +  dx) , 
given  that  there  has  been  no  failure  up  to  age   x. 

Express  the  hazard  as 


.  ,  v     dF/dx 

h(x)  =  i  -  f(x: 


h(x)dx  =  1  _dp(x)  =  "  d{log[l  -  F(x)  ]}  , 


it  then  follows  after  integration  that 
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X 

F(x)  =  1  -  exp[-  /   h(y)dy]  .  (2.3) 

0 


Thus  if  the  hazard  is  specified,  so  is  the  distribution 
function,  and  conversely. 
Note  that  if 

h(x)  =  X  >  0  ,  0  <_  x 

then 

F(x)  =  1  -  e  AX  ,       0  <_  x  ,  (2.4) 

so  a  constant  hazard  function  implies  the  exponential  dis- 
tribution of  the  random  variable   X,  and  conversely. 

Obviously  a  constant  hazard  representation  does  not 
describe  the  bath  tub  hazard  shape  of  Figure  1,  nor  does  it 
represent  a  situation  in  which  hazards  decline,  possibly 
because  design  defects  or  "bugs"  are  occasionally  removed. 
Here  are  two  hazard  representations  likely  to  be  useful  for 
such  purposes. 

1.   A  Bath  Tub  Model 

Define  the  random  variable   Z   in  terms  of  X,  X 
being  exponentially  distributed  with  mean  A   ,  as  follows: 

Z  =  G(X)  =  XL(X)R(X) 

or  (2.5) 

=  X(j)(X)  ,   cj>(X)  =  L(X)R(X) 

where 

a)  L(x)   is  concave  in  x,  L(0)  <  1,  L(°°)  =  1, 

b)  R(x)   is  convex   in   x,  R(0)  =  1,  R(0)  >  R(~) . 
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Then  the  hazard  of   Z   may  be  made  to  exhibit  a  bath  tub 
shape,  as  in  Figure  1,  by  proper  choice  of  the  functions   L 
and  R. 

Example.   Suppose 


L(x)  =  TT^        a  >  0,  0  £x 


RCx>  =  r4w     S  >  °  ' 


(2.6) 


Clearly, 


z  =  xL(x)R(x)  =  x    aX  ax 


1  +  ax   1  +  3x  1  +  ax   1  +  $x 

is  a  monotonically  increasing  function  of   x.   Furthermore, 
choose   a   large  (e.g.  a  =  10)  and   8   small  (e.g.  8  =  10   ) 
Then  it  is  intuitively  clear  that   (i)  small  x-values 
transform  into  even  smaller  z-values,  e.g.   x  =  1   corre- 
sponds to   z  =  0.91   and   x  =  2   corresponds  to   z  =  1.90, 
but   (ii)  this  effect  dwindles  as   x   increases,  so   x  =  10 
corresponds  to   z  =  9.8  and   x  =  50  to   z=47.5   and  the 
z-values  closely  resemble  the  x's  percentage-wise,  but 
(iii)  as   x  increases  still  further  the   z's  do  not  follow 
suit:   x  =  10    corresponds  to   z  =  500.   This  suggests 
that  if   x   is  a  value  assumed  by   X,  that   Z   shares  the 
properties  of   X   in  mid-range,  i.e.  for  intermediate 
x-values,  but  differs  from   X  by  having  a  disproportionate 
probability  of  assuming  small  values  (near  zero) ,  or  large 
values  (near,  but  less  than,  1/8) .   Thus  the  hazard  of 
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Z  will  appear  to  be  a  "bath  tubbed"  version  of   X, 
particularly  if   X   is  exponential. 

We  focus  attention  on  the  representation  (2.6)  in  what 
follows,  mainly  for  analytical  and  computational  convenience, 
Of  course  there  are  many  other  possibilities,  such  as 

L(x)  =  1  -  e"aX 

-6x  (2*7) 

R(x)  =  e  yx  ; 

these  latter  may  be  adjusted  to  provide  sharper-edged  tubs 
than  can  (2.6)  ,  but  iteration  of  (2.6)  may  be  induced  to 
accomplish  the  same  purpose. 

2 .   A  Decreasing  Failure  Rate  Model 

Define  the  random  variable   W   in  terms  of   X,  X 
again  being  exponential  with  parameter   ^   : 

W  =  XT(X)  (2.8) 

where   T(x)   is  an  increasing  function  of   x,  L(0)  =  1. 
Then  the  hazard  may  be  made  to  exhibit  a  decreasing 
behavior. 

Example.   Suppose 

T(x)  =l  +  cx,      c  >  0,  0  £  x  .         (2.9) 
Then 

z  =  x(l  +  ex)  (2.10) 

is   monotonic,    and   small   x-values    lead  to   comparable   z-values 


15 


(especially  when   c   is  small) ,  but  larger  x-values  are 
"amplified"  by   1  +  ex  to  yield  increasingly  large  z -values. 

Attention  will  be  focused  upon  (2.9),  although  other 
possibilities  exist  that  accomplish  the  same  purpose,  namely 
that  of  lengthening  the  right  tail  of  the  distribution  of 
X   (simulating  outliers,  for  instance)  while  leaving  the 
body  of  the  distribution  virtually  unchanged. 

B.   MATHEMATICAL  PROPERTIES  OF  THE  "BATH  TUB"  HAZARD  MODEL 

Various  analytical  properties  of  the  previously  described 
models  will  now  be  recorded.   These  provide  useful  insights 
into  the  behavior  of  the  random  variables   Z   and  the  under- 
lying (generating)  variables   X. 

1.   Monotonicity ;  Quantiles 

It  is  convenient  to  focus  on  monotonic  increasing 
transformations,  i.e.  if 

z  =  G(z)  =  Z(J>(z)  (2.10) 

then  in  order  that  the  above  function  be  monotonically  in- 
creasing, dz/dx  >  0.   Observe  that  logarithmic  differen- 
tiation of  (2.10)  provides 


d^=dx+V(x) 
z     x     ijj  (x) 


and  thus      dz/dx   >    0      if   and   only   if 


i  +  4li2E>    >  o  (2.12) 

x  4>(x) 
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Alternatively,  the  condition  is,  in  terms  of   L(x)   and 
R(x)  , 

I  +  L'  (x)  .  R'  (x)  >  Q  n    ... 

x  +  "TTxT  +  "ITUT  >   °  (2'13) 

It  is  easily  seen  that  the  important  example  (2.6) , 

jl,    i      ax        1 

9W  "  1  +  ax  '  1  +  Sx  ' 

yields  a  monotonic  relationship  between   z   and  x.   The 
fact  that  this  transformation  can  be  easily  and  explicitly 
inverted  (solved  for   x  in  terms  of   z)   will  be  exploited 
subsequently. 

Of  course  if   z(x)   is  monotonically  increasing  then  so 
is   x(z) ,  the  inverse  function.   The  events   (Z  £  z)   and 
(X   x(z))   are  equivalent,  and  so 

P{Z  <_  z]    =   P{X  <_   x(z)}  ,  (2.14) 

from  which  it  follows  that  if  x   =  x(p)   is  the  p*100% 

P 

quantile  of   X,  i.e. 

P{X  £  x(p) }  =  p  ,  (2.15) 

then 

P{Z  <_  z(p)}  =  P{Z  <_  z(x(p))}  =  p         (2.16) 

and  so   z (p) ,  the  p*100%  quantile  of   Z   is  simply  obtained 
from 

z(p)  =  x(p)  4>(x(p))  =  x(p)  L(x(p))  R(x(p))      (2.17) 

In  other  words  we  very  easily  translate  from  (points  on) 
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the  inverse  distribution  of   X  to  the  inverse  distribution 
of   Z.   Explicit  representation  of  the  distribution  of   Z 
is  however,  not  often  easily  possible. 


2 .   Hazard  and  Density  Function  Relationships 

In  order  to  investigate  the  relationship  between 
the  hazards  of   Z   and   X,  begin  by  writing 


p  =  F  (x(p)  )  =  1  -  exp[-  / 

0 


x(p) 


X 


h  (u)du] 
x 


or 


x(p) 


h  (u)du  =  -  Jtn(l-p) 


(2.18) 


Now  differentiate  with  respect  to   p   to  find 


,  ,  ,  ,,  dx(p)     1 
h  (x(p))  ,   *      =  t — 
x   ^    dp     1-j 


or 


(2.19) 


Vx(p)}  =  diV  '  Fp  =  fx(x(p})  '  ik  ;        (2*20) 


here  h    and   f   are  the  hazard  and  density  functions  of 
xx 

the  r.v.   X.   The  relationship  (2.20)  holds  for  any  distri- 
bution, of  course. 

Differentiation  of  (2.5)  reveals  the  connection  between 

h   and  h  .   From  (2.11) 
z        x 


^bL-.(p) 


1     4>'  (x(p)) 
x(p)     <Mx(p)) 


dx(p) 
dp 


=  [(j)(x(p))  +  x(p)  <J>'  (x(p))  ] 


dx(p) 
dp 


(2.21) 
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From    (4.11),    applied  now  to  the    z-hazard,    there   results 

(2.22) 


h    (z(P))    =  h    (xtp))     ^(X(P})  +X(P>    *'U(P))) 


so 


hz(z(p))    =  hx(x(p))    t<x(p))    +x(p[    V(x(p))  (2.23) 

Multiplication  of  both  sides  by   1-p  then  shows,  in  view 
of  (4.11),  that  the  density  functions  are  similarly  related 

fz(z(p))  =  fx(x(p))  ^(p))  +  x(pj-  0.(X(P))  • 

Example. 

X   is  exponential (A ) •   Then 


hz(x(p))  "  4>(x(p))  +  x(p)  cD'  (x(p)) 


Now  use  the  specific   <j)(x)   of  (2.6)  : 


ax 


*{x)    =  1  +  ax  *  1  +  3x   ' 


or,    in   terms   of    logarithms, 


In   $(x)    =    £n   ax  -    ln(l   +   ax)    -    £n(l  +   3x)       , 


so 


4>'  (x)    m   1 
<J>(x)         x 


a 


3 


1    -    agx 


(2.24) 


1   +   ax        1  +   3x        x(l   +   ax) (1   +   3x) 


,     (2.25) 


and 


<|>(x)     +    X<j>'  (x)     =    (J)(X) 


2    +    (a   +    3)x 
(1   +   ax) (1   +    3x) 


(2.26) 


finally 
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hz(z(p))  "   ax(p)  [2  +  (a  +  3)  x(p)]  (2'27) 


Although  this  expression  is  not  quite  explicit,  qualitative 
properties  of  h   can  be  deduced  from  it. 

If   p  *  0 ,  x(p)  =  -  j   £n(l-p)  i   0,  and  hence 

^'^'-w  (2-28> 


or 


lim  x(p)  h  (z(p))  =  ~-  (2.29) 

p  •*•  o 


Since  for  p  *  0 , 


and  hence 


2 

z  (p)  ~  ax  (p) 


x(p)  ~  [z(p)/a]1/2  (2.30) 

there  results 

h  (z(p)) -TTT~ — 

z         2[z(p) ]±/Z  /a 

or  , 

lim   ZzTpT  h  (z(p))  =  — - —         (2.31) 

p  ■*■  0  2/a~~ 


This  shows  that  h  (z)  t  «  as   z  ■»■  0 ,  creating  the  left- 
hand  end  of  the  bath  tub  of  Figure  1. 

If  p  ■*■   1,  x(p)  t  °°,  and   z(p)  +  1/3   so 


2   2 

v,  z,/,,!!    i  a3   x  (p) 
hz(z(P)}  ~  X   :(a  +  3) 


or 
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For   p  ■*-  1 


so 


and  thus 


l  ft2 

lira   — T  h _(z(p))  =  A        (2.32) 

P  +  ~   [x(p)]2   z  a  +  e 


1  "  *«<P>  -^HpT  (2'33) 


X(P}  -4fl  -  gz(p)  (2'34) 


h  (z(p))  -Al^-i"1)  ± 2 

(i  -  ez) 


Once  again  it  appears  that  the  hazard  rises  rapidly,  this 
time  as   x(p)  f  °°  and   z(p)  t  3   ;  the  other  end  of  the 
bath  tub  is  thus  fashioned. 

If  p  =  1  -  e   ,  then   x(p)  =  X~    .   Then 


h    (z(l  -  e-1))    -      H  +   «/X]2    [1   ±   3/A]2       (2>36) 
2  a/X    [2   +    (a   +   3)/X] 


The  bath  tub  effect  is  presumably  achieved  by  choosing   a 
large  and   3   small.   Let   a  ■*-  0   and   $  ■*•   0   independently 
in  (2.36);   it  is  clear  that  the  limiting  value  of  the 
hazard  is  X  .      This  indicates  that  the  hazard  is  (approxi- 
mately) X   for  middling  values  of   z. 

3 .   An  Explicit  Formula  for  a  Hazard 

The  expression  (2.6)  leads  to  the  relationship 

ax2(  ) 
Z(P}  =  (1  +  ax(Pn  (1  +  3x(p))  '  (2*37) 
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and  the  latter  may  be  explicitly  inverted  by  solving  a 
quadratic  equation.   The  result  is 


,_  1  _  (g  +  g)z(p)+  Aa+S)2  z2(p)+4z(p)  a(l-gz(p))   ,-  ,Q. 
(p; 2a(l  -  Sz(p)) (2'38) 


Now  a  direct  differentiation  of  this  expression  and  invo- 
cation of  (2.20)  produces  the  expression 


hz(z) 


2a(l-6z) 


6{(a+g)z+    v/(a+3)2z2    +   4ctz(l-gz) 

1   -    Sz 


+    (a   +   3) 


+ 


{  (a-B)2z  +  2g}  v/(a+6)2z2  +  4az(l-gz) 
(a   +   6)2z2    +   4az(l   -    3z) 


(2.39) 


This  form,  while  explicit,  provides  no  particularly  useful 
insights;  the  bath  tub  end  shapes  already  noted  in  (2.31) 
and  (2.35)  can  be  deduced  directly  from  (2.39). 

Some  graphical  plots  of   h    are  presented  below.  They 
illustrate  the  behavior  of  the  present  hazard  representation 
in  a  more  understandable  fashion  than  does  the  formula 
itself. 

4.   An  Explicit  Formula  for  the  Failure  Time  Distribution 
Because   x   and   z   are  monotonically  related 
through  (2.37)  we  have 
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P{Z    <    z} 


_   D  1  v   <    _  .      (a+3)z    +  Aa+$)2z2    +   4az(l-Sz)    ( 
_  P      x  <    x  2a(l   -    3z) | 

=   1   -   exp  j  -A 


+  >/T 


(a+3)z    +  v(a+3)    z      +    4az(l-3z) 
2a(l   -    3z) 


(2.40) 


Again  the  explicit  formula  seems  unproductive  of  insights. 

C.   MATHEMATICAL  PROPERTIES  OF  THE  DECREASING  FAILURE 

RATE  MODEL 

1.   The  Hazard  Behavior 

The  expression  (2.23)  can  be  applied  to  deduce  the 
hazard  function  of  the  representation  (2.8),  advertised  to 
produce  a  decreasing  failure  rate.   There  we  specified 


<|>(x)  -  T  (x)  =  1  +  ex  , 


(2.41) 


and  thus,  from  (2.2  3) 

hw(w(p))  =  (1  +  cx(p))  +  x(p)c  =  1  +  2cx(pi 

Qualitative  properties  follow  easily. 

If   p  -»■  0 ,  x(p)  4-  0,  and 

h  (w(p) )  ~  X 
w   e 

Thus  the  hazard  is  approximately   X   for  small  z, 

If   p  -»■  1,  x(p)'  t  °° ,    and 

2 
w(p)  ~  c[x(p)  ] 

so 


(2.42) 


(2.43) 
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h(w(p)) 


w 


(2.45) 


2/cw(p) 

which  clearly  decreases,  as  claimed.   It  may  be  inferred 
that  the  distribution  of  W  appears  nearly  exponential, 
but  possesses  an  extraordinarily  long  right  tail — possibly 
the  result  of  outliers. 


2.   Explicit  Formulas  for  the  Hazard  and  the 
Distribution  Function 
Direct  solution  of  the  quadratic  equation 

2 

w  =  x(l  +  ex)  =  x  +  ex 


presents 


,  .    /I  +  4cw(p)  -  1 
x(p)  =  2c~^ 


(2.46) 


which,  when  differentiated,  leads  to 


h  (w)  = 


w 


/l  +  4cw 


h..(x) 


x 


X 


/l  +  4cw 


The  distribution  function  is 

P{W£w}    =   P  |  X  £  x  =  g3ff  -    1  | 


=   1   -   exp      -X 


/I   +   4cw  -    1 
2c 


(2.47) 


(2.48) 


This  distribution  bears  a  close  family  resemblance  to  tike 
Weibull  distribution   1  -  F(w)  =  exp{-k/w},  especially  for 
large  (right  tail)  values  of  w. 
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D.   AN  ALTERNATIVE  "BATHTUB"  HAZARD  REPRESENTATION 

The  simple  parametric  model  (2.6)  leading  to  a  bathtub- 
shaped  hazard  is  by  no  means  the  only  possibility.   We  next 
describe  another  simple  approach.   It  is  that  of  defining 
a  hazard  function  having  an  appropriate  shape,  and  then 
deducing  the  corresponding  distribution  function,  and  a 
procedure  for  sampling  from  it,  rather  than  proceeding  in 
reverse  order,  as  before. 

Let  the  hazard  be  of  the  form 

h(z)  =  g(z)  +•••+  k(z)  ,  (2.49) 

where   g(z)  >  0   is  a  decreasing  function  of   z   such  that 
lim     g(z)  =0/  and   k(z)   is  an  increasing  function  of  z, 
such  that  (preferably)  k(0)  =0   and  k(°°)  =  °°.   Such  a 
function  can  yield  a  bathtubbed  hazard. 


Example . 


h(z)  =  — |-—  +  Bz  +  X  (2.50) 


A,  B,  a,  X      all  positive. 

Clearly,  (2.50)  has  a  generally  "bath  tub-like"  appear- 
ance, since 


if 


while  for 


h(z)  >  0. 


h-  (Z)  = £ — T   +  B  (2.51) 

(z  +  a) 


-  -  +  B  <  0,   then   h' (0)  <  0        (2.52) 

a 


z  >  zQ  =  |  -  a  (2.53) 
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Detailed  behavior  is  adjustable  by  choice  of  the 
parameters.   Now  the  distribution  function  of  time  to 
failure,  Z,  is  obtained  from  (2.50) : 


(    2 
P{Z  >  z}  =  exp  -  /  h(x)dx 

(   0 


■  expj"  [{tt^  +  Bx  +  A)dx] 

=   exp      -  I     A  £n(l   +   £)    +   |  z2    +    \zl 

_    /      a      .  ,      B      2.       -Az 

=    (7-TT)      exP("  I  z   }    e 


=   F^z)    F2(z)    F3(z)     , 


(2.54) 


where 


Vz)  -\<rri  A 


F0(z)    =   exp   I  -   |   z2)  (2.55) 


F3(z)    =   e 


2 


All  of  the  above  are  recognized  as  being  the  complements  of 
distribution  functions.   In  effect,  the  distribution  of   Z 
is  that  of  the  minimum  of  three  independent  random  variables 

P{Z  >  z}  =  P{XX  >  z}-P{X2  >  z}-P{X3  >  z},      (2.56) 
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X.   having  the  distribution   F.  =  1  -  P.   (i  =  1,2,3). 
This  fact  leads  directly  to  an  easy  procedure  for  simulation 
of   Z   by  simply  obtaining  the  smallest  from  among  the 
realization  of   X,  ,  X2 ,  and  X-..   The  advantage  of  the  pre- 
vious method,  based  on  (2.6)  for  instance,  is  that  only  one 
realization — that  of  an  exponential  in  that  specific  example 
— leads  to  the  realization  of   Z.   This  is  not  only  computa- 
tionally attractive,  but  seems  to  facilitate  the  application 
of  such  Monte  Carlo  variance  reduction  techniques  as  control 
and  antithetic  variables,  cf.,  Hammersley  and  Handscomb  [5]. 
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III.   OBTAINING  SPECIFIED  HAZARD  BEHAVIOR 
BY  SIMPLE  SAMPLING 

The  development  of  the  last  section  illustrates  one 
manner  in  which  hazard  behavior  may  be  conveniently  repre- 
sented and  simulated.   We  now  show  how  such  behavior  may 
alternatively  be  obtained  by  simple  simulation,  i.e.  from 
one  realization  of  a  basic  (possibly  exponential)  random 
variable . 

Refer  to  (2.5),  in  which 

Z  =  G(X)  (3.1) 

and,  if   G(«)   is  monotonically  increasing, 

z(p)  =  G(x(p))  ,  m  (3.2) 

z(p)       and      x(p)      being   the   p-100%    percentiles   of      Z      and      X, 
respectively.      Then   the   counterpart   to    (2.2  3)    that "results 
from  differentiation   of    (3.2)    is    the   expression 

h,(«(p)J    =  hx(x(p))    G,(xip))    =  hx(x(p))    TH^gxT    (3.3) 

Consequently,  if  one  specifies   h  (z)   as  a  suitable  func- 
tion  of  the  "time"   z,  and  specifies  the  distribution  of 
the  stochastic  variable   X — and  hence  its  hazard,  h — there 
results  a  differential  equation  for   z(x)  =    G(x) : 

hz(z»  I  =VX)  >  (3-4) 

integration  then  provides  the  desired  transformation,  G.   In 
other  words,  we  seek   z(x)   satisfying 
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z  X 


/  h  (u)du  =  /  h  (v)dv  ,  (3.5) 


0    2         6 


which  can  sometimes  be  carried  out  in  a  useful  closed  form. 


Example  3.1.   Refer  to  the  example  of  Section  II,  wherein 
h   is  given  by  the  expression  (2.50)  and  we  assume  that  X 
is  exponential,  so  h   is  constant.   Then  in  order  to  deter- 
mine  G(x)  =    z(x),  solve  the  equation 


/   I  — x —  +  Bu  +  A 
q   I  u  +  a 


x 

du  =  /  dv 
0 


or 


A  fcn(l  +  -)  +  |  z2  +  Az  =   x  (3.6) 

a     2 


Closed-form  solution  of  this  expression  for   z   in  terms  of 
x   is  of  course  impossible.   One  possible  approach  is  purely 
numerical:   find  an  approximate  solution,  z~ (x) ,  e.g.  the 
appropriate  solution  of  the  quadratic 

|  z2  +  Az  =  x  (3.7) 

and  then  correct  the  result  by  a  few  Newton- Raphson  itera- 
tions.  In  other  words,  put 


,  »    -  A  +  vA   +  2Bx  /.,  Qs 

z1(x)  =  g ;         (3.8) 


now   apply  Newton    to   obtain   an   improved  solution 

A   £,n(l  +  Zj/a) 
22(X)     *    21    "   A/(z1  +  a)+Bz1+A  (3'9) 
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which  will  be  feasible  if   0  <  z2.    The  process  can  be 
iterated  (the  numerator  will  change  after  the  first  iter- 
ation) .   If  one  wishes  to  use  this  model  it  may  actually  be 
desirable  to  start  by  solving 


|  z2  +  (A  +  hz    -  x  =  0  (3.10) 

for   z.  ,  in  which  case  the  numerator  will  not  be  as  shown  in 
(3.9) ;  convergence  may  be  more  rapid. 

Example  3.2.   Change  the  hazard  representation  of  the  previous 
example  as  follows:   let 

h  (u)  = 7  +   Bu  +  A;     (A,  a,  B,  A  >  0)     (3.11) 

(v  +  a) 

then 

/   h  (u)du  =   ,  A2i r  +  |  z2  +  Az   (3.12) 

i         z         a(z  +  a)    2 

Now  it  is  necessary  to  solve 


Az    ,  B   2  ,  ,  . _  n  _ , 

+  7T   z   +  Az  =  x  ,  (3.13) 


ct(z  +  a)    2 
i.e. ,  the  cubic 

|  z3  +  [a  |  +  A]z2  +  [-  +  aA  -  x] z  -  ax  =  0  ,     (3.14) 

which  can  be  carried  out,  at  least  formally,  in  closed  form. 

Once  again  an  iterative  solution  that  begins  by  dropping  the 

cubic  term,  solving  the  resulting  quadratic  for   z, (x) ,  and 

then  continuing  along  the  Newton-Raphson  road  may  be^ 

successful.   Further  investigations  of  these  ideas  should 

be  conducted. 
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IV.   COMPUTER  SIMULATION  AND  ESTIMATION  PROCEDURE 

A.   SIMULATION  AND  NUMERICAL  RESULTS 

In  previous  chapters  an  analytical  model  was  described 
for  the  failure  rate  function;  useful  formulas  were  also 
derived  from  the  model  (2.39)  and  (2.40).   Before  the  model 
is  used  in  realistic  situations,  it  will  be  convenient  to 
build  a  computer  simulation  model  for  model  validation. 

1.  Algorithm 

First  a  very  basic  simulation  model  was  built  for 
determining  the  general  shape  of  the  failure  rate  function 
associated  with  parameters   a,  3  and   A.   In  the  simulation 
model,  a  was  selected  to  be  1.0  and  2.0,  3   was  selected  to 
be  0.05,  0.01,  0.005,  0.001  and   A   was  selected  to  be  0.1. 
These  values  were  picked  arbitrarily,  it  is  only  stipulated 
that   a  is  always  greater  than   3.   The  general  algorithm 
of  the  simulation  model  is  shown  in  Fig.  2. 

The  hazard  function  is  calculated  according  to  the 
model  (2.37)  and  the  system  logic  function  (2.40).   The 
results  were  shown  in  Fig.  3  and  Fig.  4. 

2 .  Some  Comments 

These  simulation  results  show  that: 
a.   Parameter   a   is  effective  when   z   values  relatively 
have  small  values.   That  is,  it  influences  the  early 
failure  period. 
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INITIAL  VALUES  OF 
PARAMETERS 


GENERATING  RANDOM  NUM 
FROM  EXPONENTIAL 
DISTRIBUTION 


GETTING  Z  VALUES 
Z=G(X) 


GETTING  FAILURE  RATE 
FUNCTION 


PLOTTING  Z  vs  FAILURE 
RATE  FUNCTION 


Fig.  2.   Basic  Algorithm, 
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b.  Parameter   3   is  effective  on  the  relatively  bigger   z 
values.   That  is,  it  describes  wearout  failures. 

c.  Parameter   A   has  little  effect  on  the  shape  of  the 
curve;  it  is  a  scale  factor. 

d.  The  last  important  observation  is  that  the   z   values 
are  limited  by  the  parameter   3   such  that : 

-ti 

When   z   equals  1/3/  h  (z)   goes  to  infinity. 

B.   EMPIRICAL  DENSITY  FUNCTION  AND  ESTIMATION  PROCEDURE 

In  this  section,  an  estimation  procedure  for  parameters 
a,   3/  A  is  defined  and  the  procedure  of  finding  the 
empirical  density  function  is  described. 

1.   Empirical  Density  Function 

Suppose  that   N. ,  N,  A   and  N~   are  defined  as 
follows : 


N.  =  the  frequency  of  data  points  for  each  time  interval 

between   a.   and  b.   for   i   equals  l,2,3,...,k. 
N  =  the  total  number  of  failures  where, 

k 

n  =    y    n. 
i=i    x 

A   =  the  length  of  each  interval 


A  =  b.  -  a. 
l    l 
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Nq  =  the  total  number  of  survivors  by  the  beginning  of 
each  interval  which  can  be  shown  as  follows: 

N  if   i  =  1 

N  -  I      N.         if   i=2,3,...,k 
j-1   D 


where   > .  ,  N .   is  the  number  of  failures  before 

interval   i,  or  N  -  V.  n  N.   is  the  number  of 

^=1  j 

failures  after  interval   i. 


If   h  (z)   is  defined  as  the  density  function  of 
z  ■* 

hazard,  then  the  relationship  between  the  frequences  of  the 
failure  data  and  density  function  of  hazard  can  be  given 

as  follows: 

b. 

l 

**.;£  NQ  •   /   hz(z)dz  (4.1) 

a. 

i 

where   i  =  l,2,3,...,k   and   a.  =  (i-l)A,  b.  =  iA.   At  this 
point  some  approximation  can  be  made  in  expression  (4.1), 

such  that 

b. 

l 

/   h  (z)dz  ~  (b.  -  a.)-h  (i(b.  -  a.))  (4.2) 

*  z       ~    1     l    z     1     1 

ai 

for  small  values  of   A  =  b-a.   Then  this  approximation  (4.2) 
is  substituted  in  the  expression  (4.1),  it  turns  out  as 
follows : 


or 


N.  =  A-N.  h  (iA)  (4.3) 

l       0   z 
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N. 
hz(iA)    =   N~~^A       '  i   =   1>2"    "k  (4*4) 

for   the   general   case,    it  will  be: 


N. 

1 

N 


for      i  =   1 


h2(iA)    =    ^  (4.5) 

Ni  1 


i-1  -A 

(N  -     I      N.)  for      i   =    2,3,  .  .  . 


j-l 
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The  expression  (4.5)  will  be  used  in  the  calculation  of 
an  empirical  hazard  function  and  also  used  in  the  proposed 
estimation  procedure. 

2 .   Estimation  Procedure  for  Model  Parameters 

Two  approaches  can  be  used  for  this  problem.   The 
first  idea  is  to  use  the  relationship  between  a  sample  pth 
percentile  and  the  related  probability  of  the  pth  percentile. 
A  subsequent  idea  is  to  approach  the  problem  as  a  nonlinear 
least  square  estimation  problem  for   a,   3a   A:   pick   a,  3 
A   so  that  the  values  obtained  minimize  the  sum  of  squared 
errors  in  an  objective  function. 

a.   3-Percentile  Approach 

Suppose   Fv(x)   is  the  cumulative  probability 

A 

function  of  the  exponential  distribution.   The  pth  percentile 
x(p)   is  equal  to  the  value  such  that: 

p  =  Fx(x(p))  =  1  -  e'Xx(p)  (4.6) 
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Then 

x(p)  =  F_1(p)  >  -  i   £n(l-p)  (4.7) 

where   A  >  0 . 

If   z(p)   is  defined  to  be  the  pth  percentile 
in  the  bathtub  model  (2.37)  then   z(p)   and  x(p)   have  a 
definite  relationship  with  expression  such  that: 


-<n\  =  a*  (p) 

^P;    (1  +  ax(p))  (1  +  3x(p)) 


by  substituting  (4.7), 


a{-  i  £n(l-p) }2 


z(p)    = 


1  + 


a 


j  An(l-p) 


1  +  3 


or 


-  -  £n(l-pi 


(4.8) 


z(p)    = 


A2)  £2( 


I1  +  Xe<P>}{  x  +X£(P}} 


(4.9) 


where   e (p)  =  -  Zn(l-p). 

In  the  last  equation  (4.9)  ,  there  are  three  un- 
known variables  which  are   a,   3,   A.   If  three  independent 
equations  associated  with  expression  (4.9)  are  defined, 
clearly,  they  can  be  solved  for  the  unknown  parameters. 
Actually,  this  can  be  easily  done,  for  using  different  values 
of  percentiles  such  that: 
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(a/A2)s2(P1) 
z(P,)    = 


1      {1+fe(?i)}{1  +  xe(pi) 


(a/A2)e2(p.) 
z(P2)  = 


{l+f  £(P2)}{l  +  |£(P2)} 


(a/A2)e2(p,) 
z(P3)  = 


{l   ^MP3)}{ltfe(P3)} 


where   z(P.)   and   e(P.)   are  the  known  values  associated 
with   P..   Actually,  to  get  the  percentile  values,  two  kinds 
of  approaches  can  be  made.   First  they  can  be  computed 
directly  from  data  using  the  simple  statistics  method,  such 
as  z/-\       is  (approximately)  z(i/(n+l)).   Second  they  can 
be  computed  from  the  empirical  density  function.   Generally 
the  choice  depends  on  the  form  of  the  data  that  are  avail- 
able. 

To  decide  for  the  effectiveness  of  this  type  of 
estimation,  another  basic  computer  simulation  is  made  by 
modifying  the  first  computer  simulation  algorithm.   This 
algorithm  is  shown  in  Fig.  5. 

The  simulation  was  run  four  hundred  times  and 
in  each  run  the  sample  size  was  assumed  to  be   n  =  50. 
Initially  the  parameters   a,  8   and  X  were  taken  to  be, 
respectively,  1.0,  0.05  and  0.1.   The  percentiles   P,,  P2 
and   P.   were  used  in  each  run  such  that  0.1,  0.5  and  0.9 
respectively.   The  estimation  results  are  shown  in  Table  I. 
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SELECTION  OF  DIFFRENT 
VALUES  OF  PARAMETERS 


GENERATING  RANDOM  NUM 
FROM  EXPONENTIAL 
DISTRIBUTION 


GET  THE  FREQUENCY 
DISTRIBUTION  FOR 
CERTAIN  INTERVAL 


FIND  EMPIRICAL  DENSITY 
FOR  HAZARD, ESTIMATE 
THE  PAREMETERS 


Fig . 5 .Simulation  algoritm  for  estimation  procedure 
and  emprical  density  function. 
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As  can  be  seen,  the  estimated  parameters   a,   3, 
X        have  quite  different  values  in  all  runs.   To  judge  their 
sampling  distributions,   histograms  were  drawn  separately; 
Figures  6,  7  and  8  are  the  histograms  of  the  estimated 
values  of   a,   3/    A,  respectively. 

It  is  suggested  by  Figures  6  and  7  for  A   and   3, 
that  normality  may  be  assumed  with  sample  means  0.9246, 
0.05488  and  sample  variances  0.001065,  0.0002246,  respec- 
tively.  Their  mean,  median,  trimean  and  midmean  are  pretty 
close  to  each  other.   But,  in  the  case  of   a,  the  histogram 
(Fig.  8)  is  a  quite  different  picture  which  looks  something 
like  an  exponential  distribution  instead  of  normal  distri- 
bution.  The  reason  is  clear,  because  the  negative  values 
of  estimated   a  were  not  taken  into  account.   They  are 
physically  infeasible. 

b.   Nonlinear  Least  Squares  Approach 

The  least  squares  criterion  used  here  can  be 
stated  formally  and  generally  as  follows: 

N        A   2 
minimize    £   (Y.  -  Y.)  (4.10) 

i=l    1    X 

where   N   is  the  number  of  observations,  Y.   is  the  fitted 
value  of   Y..   In  this  case,  expression  (4.10)  can  be  re- 
written as  follows 

min  I    (Z.  -  Z. (a,3,A)}z  (4.11) 

i=l   1 
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TABLE  I 
COMPUTER  SIMULATION  RESULTS  FOR  ESTIMATION  PROCEDURE 


Sample  size   n  =  50; 
Value  of  parameters  a,  3, 


Number  of  Total  Runs  m=400; 
X  Introduced:  a=l,  3=0.05,  A=0 . 1 


Number  of 

Runs 

A 

a 

3 

1 

0.0858 

2.2255 

0.0549 

2 

0.06591 

1.2745 

0.0779 

3 

0.08381 

9.498 

0.0688 

4 

0.0881 

0.9304 

0.039 

5 

0.12  86 

3.7429 

0.0498 

6 

0.0883 

2.114 

0.069 

7 

0.0785 

2.8281 

0.0593 

8 

0.0945 

0.4073 

0.0584 

9 

0.089 

0.3833 

0.0452 

10 

0.1344 

1.3351 

0.0426 

11 

0.0677 

0.2111 

0.0  729 

12 

0.1304 

2.4871 

0.0512 

13 

0.0714 

0.3211 

0.0623 

14 

0.069 

0.3547 

0.0657 

15 

0.0775 

0.4813 

0.0772 

16 

0.1618 

2.7416 

0.027 

17 

0.0827 

0.927 

0.0705 

13 

0.0816 

1.1663 

0.0836 

19 

0.0515 

0.3777 

0.0775 

20 

0.10  89 

1.6439 

0.0522 

21 

0.0767 

0.8667 

0.0605 

22 

0.0499 

0.3414 

0.0748 

23 

0.0569 

0.2118 

0.0608 

24 

0.0648 

0.2689 

0.0715 

25 

0.1415 

2.1178 

0.0518 

26 

0.16  36 

1.3635 

0.0186 

27 

0.0947 

6.101 

0.0664 

42 


TABLE    I    Cont. 


Number  of 
Runs 

A 

a 

8 

28 

0.1063 

0.5235 

0.0623 

29 

0.0825 

2.6602 

0.059 

30 

0.13 

0.3412 

0.0303 

31 

0.817 

0.24 

0.0685 

32 

0.0843 

0.4377 

0.0393 

33 

0.1445 

4.6941 

0.0442 

34 

0.0737 

0.2956 

0.0665 

35 

0.0  45 

0.1372 

0.0646 

36 

0.1108 

1.1026 

0.0427 

37 

0.1277 

1.091 

0.052 

38 

0.0514 

0.1362 

0.062 

39 

0.1165 

0.6  79  4 

0.0469 

40 

0.0954 

3.7845 

0.0587 

41 

0.0563 

0.4123 

0.0738 

42 

0.0684 

0.2903 

0.0574 

43 

0.1263 

2.059 

0.0357 

44 

0.1335 

7.1595 

0.0226 

45 

0.0708 

0.465 

0.063 

46 

0.0872 

2.16  4 

0.0657 

47 

0.0834 

0.3568 

0.0567 

48 

0.0429 

0.127 

0.0694 

49 

0.0428 

0.1295 

0.0639 

50 

0.0834 

0.6807 

0.0489 

51 

0.1141 

1.5349 

0.0566 

52 

0.1057 

1.1757 

0.0441 

53 

0.110  4 

5.0692 

0.0433 

54 

0.1426 

0.8845 

0.0244 

55 

0.10  84 

0.5775 

0.0505 

56 

0.1048 

1.8399 

0.0562 

57 

0.0561 

0.3197 

0.0643 

58 

0.0607 

0.4564 

0.0607 

59 

0.0724 

0.3133 

0.0551 
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TABLE    I    Cont. 


Number   of 

Runs 

A 

a 

6 

60 

0.0957 

0.9584 

0.0557 

61 

0.0804 

0.7378 

0.0423 

62 

0.0888 

0.1718 

0.0393 

63 

0.064 

0.3779 

0.9543 

64 

0.0324 

0.2202 

0.0791 

65 

0.0695 

0.2428 

0.0564 

66 

0.0903 

6.4194 

0.0542 

67 

0.0741 

0.4118 

0.0718 

68 

0.0462 

0.167 

0.0671 

69 

0.0805 

1.0279 

0.0662 

70 

0  .0  469 

0.2157 

0.06.9  3 

71 

0.0487 

0.2457 

0.0688 

72 

0.0682 

1.1847 

0.0665 

73 

0.0538 

0.1964 

0.0573 

74 

0.1516 

2.099 

0.0352 

75 

0.0  76  3 

0.3067 

0.0679 

76 

0.0558 

1.5589 

0.0  79 

77 

0.0463 

0.2252 

0.0612 

78 

0.095 

5.6243 

0.0627 

79 

0.190  3 

6.0 

0.0604 

80 

0.0822 

0.9989 

0.0738 

81 

0.066 

0.3422 

0.0926 

82 

0.0784 

0.3336 

0.0517 

83 

0.0499 

0.1453 

0.062 

84 

0.0667 

0.1357 

0.064 

85 

0.0431 

0.0761 

0.0568 

86 

0.1485 

6.1674 

0.0262 

87 

0.115  4 

1.0451 

0.0476 

88 

0.0043 

0.3493 

0    0648 

89 

0.0983 

1.3101 

0.0422 

90 

0.138 

6.7408 

0.0443 
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TABLE    I    Cont 


Number   of 
Runs 

X 

a 

0 

91 

0.1117 

0.7281 

0.058 

92 

0.065 

0.4284 

0.0482 

93 

0.0638 

1.2754 

0.0611 

94 

0.1603 

1.003 

0.0302 

95 

0.0436 

0.1647 

0.085 

96 

0.0905 

6.2455 

0.0  70  3 

97 

0.1309 

1.4348 

0.0431 

98 

0.0787 

0.5245 

0.0686 

99 

0.1927 

7.801 

0.009 

100 

0.0683 

0.525 

0.0698 

101 

0.1163 

1.7146 

0.0355 

102 

0.1005 

1.5  79  5 

0.0463 

103 

0.1492 

1.5169 

0.0465 

10  4 

0.1511 

2.7243 

0.318 

105 

0.0  82  3 

0.1634 

0.0438 

106 

0.1493 

6.9467 

0.0437 

107 

0.1185 

9 .4492 

0.0603 

108 

0.10  39 

0.589 

0.0406 

109 

0.0832 

0.3935 

0.0612 

110 

0.059 

0.145 

0.0508 

111 

0.133 

8.6071 

0.0  39  7 

112 

0.0877 

0.4751 

0.0448 

113 

0.0772 

0.2752 

0.0624 

114 

0.0659 

0.4782 

0.0583 

115 

0.0796 

0.2978 

0.0537 

116 

0.112  8 

1.668 

0.0601 

117 

0.0974 

1.2931 

0.0515 

118 

0.1105 

0.3394 

0.0468 

119 

0.1216 

1.0529 

0.0655 

120 

0.1046 

2.0384 

0.0591 

121 

0.128 

i ■ 

3.2468 

0.0405 
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TABLE    I    Cont. 


Number   of 

Runs 

A 

a 

6 

122 

0.0818 

1.092 

0.0584 

12  3 

0.0924 

2.4295 

0.0439 

124 

0.0732 

0.1916 

0.0495 

125 

0.1284 

1.8515 

0.0407 

126 

0.1376 

2.6651 

0.0282 

12  7 

0.0877 

7.8429 

0.0703 

12  8 

0.0865 

1.1942 

0.0438 

129 

0.1404 

2.8991 

0.0326 

130 

0.0452 

0.0978 

0.0663 

131 

0.0928 

3.046  3 

0.0471 

132 

0.0668 

0.4602 

0.0551 

133 

0.2089 

0.7373 

0.0119 

134 

0.0946 

0.4011 

0.038 

135 

0.1352 

1.2162 

0.0284 

136 

0.0597 

0.2241 

0.0539 

137 

0.10  6 

1.7821 

0.0  35  4 

138 

0.0508 

0.2689 

0.0596 

139 

0.0522 

0.1244 

0.055 

140 

0.10  4 

0.8464 

0.0603 

141 

0.1002 

0.817 

0.0711 

142 

0.0692 

1.4214 

0.0744 

143 

0.0995 

0.5442 

0.0435 

144 

0.0697 

0.4966 

0.0487 

145 

0.0737 

0.8995 

0.0708 

146 

0.059 

0.2998 

0.0654 

147 

0.1235 

1.1251 

0.045 

148 

0.1023 

0.615 

0.0349 

149 

0.1 

0.9087 

0.0466 

150 

0.1242 

2.9267 

0.0356 

151 

0.0966 

0.4989 

0.0616 

152 

0.0913 

2.5416 

0.0638 

15  3 

0.1301 

1.3897 

0.0434 

46 


TABLE  I  Cont. 


Number  of 
Runs 


a 


3 


154 
155 

15  6 
157 
158 
159 
160 
161 
162 

16  3 
164 
165 
166 
16  7 
168 
169 
170 
171 
172 
173 
174 
175 
176 
177 
178 
179 
180 
181 
182 
183 
184 


0.1033 

0.0638 

0.1408 

0.1311 

0.0608 

0.0574 

0.0998 

0.1114 

0.0415 

0.1586 

0.0723 

0.0984 

0.1011 

0.093 

0.1141 

0.1953 

0.0568 

0.0927 

0.0632 

0.0536 

0.1033 

0.0633 

0.1414 

0.10  89 

0.1165 

0.1731 

0.0679 

0.1682 

0.0868 

0.0  80  7 

0.0778 


2.0535 

0.2219 

3.6331 

4.9526 

0.469 

0.5886 

4.912 

0.4898 

0.1827 

7.5257 

0.9055 

2.1521 

2.1535 

4.8488 

4.2636 

4.6512 

0.8  36  4 

5.1303 

0.7399 

0.2732 

0.9515 

1.6597 

1.2547 

4.1786 

0.5441 

4.2913 

0.1064 

3.1417 

0.6399 

1.046 

1.185 


0.0466 

0.0594 

0.0434 

0.0495 

0.0462 

0.0723 

0.0554 

0.0363 

0.0854 

0.0397 

0.0637 

0.0675 

0.0573 

0.0555 

0.0443 

0.0452 

0.07099 

0.068 

0.069 

0.07138 

0.0787 

0.0791 

0.0318 

0.051 

0.0432 

0.0192 

0.0451 

0.0332 

0.0685 

0.0607 

0.0694 
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TABLE    I    Cont 


Number  of 

Runs 

X 

a 

6 

185 

0.0769 

0.1999 

0.0009 

186 

0.0396 

0.119  4 

0.064 

187 

0.1012 

0.499 

0.485 

188 

0.1389 

9  .  70  85 

0.0334 

189 

0.0665 

1.1661 

0.0599 

190 

0.0471 

0.1374 

0.0643 

191 

0.1316 

4.1676 

0.0259 

19  2 

0.0687 

0.6692 

0.0928 

19  3 

0.0866 

0.5419 

0.0484 

194 

0.0966 

1.579 

0.0524 

195 

0.0773 

1.0102 

0.0668 

196 

0.1076 

0.7413 

0.0704 

197 

0.0  85  7 

0.8715 

0.0548 

198 

0.0  80  3 

1.110 

0.0549 

199 

0.052 

0.1396 

0.0598 

200 

0.0757 

0.9558 

0.0622 

1 
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The  goal  is  to  choose  values   a,  8   and  X  which  minimize 
the  expression  (4.11) .   In  the  linear  case,  this  can  be  done 
using  elementary  calculus  by  taking  the  partial  derivatives 
with  respect  to   a,  3   and  X    ,    setting  each  of  them  equal 
to  zero  and  solving  the  resulting  linear  equations  (Normal 
equations) .   But  the  present  problem  is  different  because 
the  expression  (4.11)  is  highly  nonlinear  in  the  parameters, 
and  indefinite  in  terms  of  convexity.   This  is  why  a  non- 
linear programming  technique  is  used. 

There  are  a  few  general  approaches  to  the  solu- 
tion of  the  nonlinear  estimation  problem.   One  of  them  is 
the  direct  optimization  approach. 

Specifically,  equation  (4.9)  is  considered  as 


follows : 


PXP2£:     (p) 

z(p)   =  11  +  p,s(P)Hl  +  fT^TpTI  (4*12) 


where 


pl   =  I 


p2    =   X 


p3    =   I 


In  equation    (4.12),    z(p)       and      e  (p)       are  known  values,    p^^ 
are    unknown   variables.      If   it   is    rewritten   as    the   sum  of 
squared  errors   the   objective    function   will   be: 


j  P]_P2£    (Pi) 


Xx      |z(Pi}    "    (1   +    p2e(P.))(l   +    p3£(Pi))     (  (4*13) 
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Now  nonlinear  optimization  problems  can  be  stated  such  that 

N  (  PlP2£2(P  ) 

^n  Jx  j  Z(Pi}  *  (1  +P2e(P.))(l  T  p3£(P.)) 

Subject  to  p,  >_  0 ,      p2  >  0,   p3  >_  0 

In  order  to  perform  the  optimization,  i.e.  to 
solve  this  nonlinear  optimization  problem,  the  GRG  package 
(Generalized  Reduced  Gradient)  [6]  was  used;  it  is  very 
convenient  for  this  case.   Ten  runs  were  made  with  five 
different  initial  points,  both  with  analytically  computed 
derivatives  and  with  numerically  computed  derivatives.   The 
first  four  initial  points  were  selected  arbitrarily.   The 
last  one  was  picked  such  that: 


Pn  =  ~ 

1  A 

P0  =  -4^~                                                      (4-14) 

2  A 


p3  = 


Z(50) 


where   zm   is  tne  smallest  number  in  the  sample,  z  i^q) 
is  the  largest  number  in  the  sample,  and   A   is  the  natural 
logarithm  of  p  =  0.5   divided  by  the  median  of  the  sample. 
Initial  points  are  shown  in  Table  II. 

The  optimization  results  of   a,  3,   and   p 
associated  with   p.,,  p2   and   p3    are  tabulated  in 
Table  III  and  Table  IV,  respectively. 
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As  it  is  seen  from  Tables  III  and  IV,  the  program 
was  terminated  at  different  optimal  points.   Because  the 
objective  function  is  highly  nonlinear  and  apparently  in- 
definite (neither  convex  nor  concave) ,  it  may  sometimes 
have  stopped  at  a  local  minimum  rather  than  at  a  global 
optimum  point.   But  if  the  results  are  compared  to  the 
3-percentile  results,  the  nonlinear  least  square  estimates 
show  much  greater  accuracy  and  seem  more  consistent. 


TABLE  II 

INITIAL  GUESS  POINTS 

#  of 

Points 

pl 

p2 

p3 

I 

50.0 

2.0 

30.0 

II 

10.0 

10.0 

0.5 

III 

80.0 

80.0 

5.0 

IV 
V 

1.0 

1.0 

1.0 

According  to  the  e< 

3.  (4.14) 

54 


CO 


<C 

D 

a 

CO 

E-" 

CO 

< 

w 

J 

tf 

4 

w 

2 

H 

kI 

2 

O 

2 

U 

in 

2 

• 

H 

o 

CO 

D 

il 

>H 

co 

H 

CQ 

Q. 

H 

H 

n 

N 

w 

a 

O 

J 

o 

<3 

a 

^H 

Eh 

2 

< 

II 
CN 

CN 

a. 

Q. 

>k 

.» 

rH 

o 

Q. 

• 

o 

fci 

rH 

O 

II 

CO 

W 

H 

D 

Q. 

3 

> 

CO 

ID 

Q 

3 

W 

rH 

En 

(tJ 

aS 

> 

H 

rH 

Eh 

(d 

CO 

a; 

W 

OS 

0) 

3 

rH 
(« 
> 

> 

•H 
-P 
O 
d) 
•co 
.Q 
O 


co 

a 


CN 

Q. 


cq 


cq 


CQ 


CQ 


rH 

rd  p 

•H  C 
P  -H 
•H  0 

c  ft 


4-1  CO 
0  § 


CT» 

00 

in 

co 

CM 

rH 

<y> 

00 

o 

O 

o 

cr» 

cr» 

o 

o 

rH 

o 

o 

rH 

rH 

in  in  in  in  in 


CO 

rH 

m 

» r* 

VD 

o 

•^r 

•^ 

■^r 

>£> 

CO 

•^r 

rH 

*J« 

CN 

00 

<N 

CN 

O 

CN1 

o 

in 

O 

rH 

r* 

T 

r^ 

rH 

rH 

in 

vo 

rH 

r*» 

rH 

V0 

CN 

m  m  m  in  in 


•^ 

r^ 

00 

10 

"=r 

CN 

CM 

CN 

CN 

-3" 

CO 

CO 

CO 

CO 

co 

CO 

en 

CO 

CO 

co 

CN   "* 


o 

CN 

in 


CN 


m 
en 


in  oo  cn 


co  in 
co  ko 


co 


r-  r-» 


V0 

m 
oo 

CN 


O 
CO 

o 


in  co 


^o  r»  ko  r-  r» 

r*-  r-»  r^  r»  r» 

r-~  r»  r^  r-»  r-» 

rH  rH  rH  rH  rH 


CO   CO   CO   CO   CO 


co  co  co  co 


o  o 


o  o 


CO   rH   CN  \£) 

o  o  o  cr» 


m  in 


in 

rH 


m 


r*  rH 


CN 


KO  10  CO 

CO  00  00 

o  o  o 

CO  CO  CO 


V0 

in 

o 
co 


in 

CN 

CO 


00   CO   CO   00   00 


co 


co 
co 

CO 


in 

in 


rH 

co 

V0 


CN 
O 


CO 


00   CO   CO   CO   CO 


H   H 
H 


>   > 


00 

o 

00 

in 


CN 

CO 

in 


rH 

oo 
in 


oo 
o 
oo 
in 


co 
in 
in 


"^  CO  00 

<J\  CXl  ^0 

-ST  CN  rH 

in  cn  in 


co 

CT> 


m  vo  in  c?>  in 

rH  rH  rH  rH  rH 

Cl  CT\  CT\  0">  Cft 

in  in  in  in  in 


CO 

<T\ 

r^ 

•<3* 

<Tt 

CN 

V0 

r- 

rH 

10 

co 

tr 

o 

V0 

CN 

V0 

00 

■H 

o 

in 

CO 

CO 

in 

in 

co 

00 

00 

o^ 

10 

<x> 

CO 

co 

CO 

CO 

CO 

■^ 

in 

•^ 

<o 

rH 

H       H       H 
H      H 


CN 


CN 

O 


m 


oo    cri    m     oo 

O       <T»       rH       o 


CO 


rH 

oo 
in 


co 


CN       CTi 


>       > 


• 

a) 

<u 

> 

> 

•H 

•H 

p 

P 

fO 

(tJ 

> 

> 

•H 

•H 

rH 

SH 

0 

CD 

T3 

T3 

tn 

P 

c 

3 

•H 

0 

W 

.c 

3 

p 

•H 

>i 

5 

-Q 

73 

T3 

<D 

<U 

P 

P 

3 

3 

a 

a 

e 

e 

0 

0 

u 

u 

<:   cq 


55 


'S' 

VD 

r^ 

rH 

ro 

CO 

ro 

CN 

r^ 

■^ 

00 

00 

CO 

CO 

00 

0) 

p- 

en 

en 

rH 

CO 

00 

•^r 

rH 

r^ 

rH 

rH 

rH 

f-i 

rH 

rH 

3 

CQ 

•sr 

*3« 

■^ 

in 

^ 

CO 

en 

en 

00 

en 

CN 

CN 

CN 

CN 

CN 

r-H 

(0 

> 

m 

CO 

ro 

ro 

ro 

•^ 

•^ 

•^ 

■^ 

•^r 

00 

00 

CO 

00 

00 

in 

LO 

in 

^n 

in 

r- 

r-» 

r^ 

r^ 

r>- 

ro 

ro 

ro 

ro 

ro 

_"■  " 

> 

•H 

m 

<£> 

0 

0 

in 

■^ 

r^ 

ro 

en 

CO 

■"^ 

CN 

ro 

0 

en 

-M 

<3« 

ro 

vo 

rH 

CN 

"?r 

r- 

00 

CN 

en 

00 

rH 

r-» 

r-« 

rH 

U 

<: 

<^> 

<Tl 

00 

VD 

CN 

en 

en 

en 

CT> 

en 

CN 

•** 

in 

0 

CN 

<U 

H 

ro 

ro 

ro 

in 

^r 

^r 

•^ 

^ 

^ 

en 

CO 

vo 

en 

CO 

■r—i 

A 

U3 

LD 

in 

in 

en 

r- 

r^ 

r-« 

r-» 

r-» 

<* 

ro 

•^ 

ro 

ro 

O 

rH 

CN 

LO 

•^r 

0 

-sT 

cn 

0 

■^ 

CN 

in 

0 

0 

0 

0 

0 

rH 

iH 

0 

CN 

<-i 

0 

0 

rH 

O 

00 

m 

rH 

vo 

vo 

00 

CQ 

PI 

ro 

ro 

ro 

m 

r>~ 

in 

vo 

00 

in 

CN 

CN 

CN 

CN 

CN 

CN 

CN 

CN 

CN 

CN 

"* 

^r 

t 

<* 

^ 

in 

in 

in 

m 

in 

O 

O 

0 

0 

0 

0 

0 

0 

O 

0 

rH 

rH 

rH 

t-^ 

rH 

ro 

Q. 

00 

•^r 

ro 

en 

<£> 

r* 

<3< 

0 

O 

0 

en 

rH 

m 

ro 

00 

^0 

CN 

CO 

0 

ro 

00 

rH 

rH 

ro 

"^ 

^ 

O 

CN 

'O 

00 

< 

^r 

ro 

CN 

ro 

r-» 

^ 

^r 

"3" 

m 

ro 

^ 

CN 

00 

CN 

CN 

CN 

CN 

CN 

CN 

rH 

^r 

•^r 

•^ 

^ 

■^ 

0 

-tf 

0 

ro 

in 

0 

0 

O 

0 

rH 

0 

0 

0 

0 

0 

rH 

rH 

rH 

rH 

rH 

r~ 

ro 

en 

■^r 

0 

rH 

r- 

O 

rH 

r- 

en 

^0 

in 

in 

^0 

CN 

ro 

CN 

CN 

ro 

O 

0 

O 

O 

en 

• 

CN 

r^ 

"3" 

00 

VD 

<H 

rH 

rH 

rH 

rH 

O 

0 

O 

O 

en 

4J 

00 

en 

0 

rH 

ro 

rH 

rH 

rH 

r-i 

c 

CQ 

O 

m 

CN 

CN 

ro 

0 

u 

<-* 

0 

rH 

0 

rH 

00 

rH 

HI 

CN 

Q. 

H 

H 

ro 

en 

*0 

O 

p* 

in 

CO 

00 

ro 

en 

O 

CO 

en 

r^ 

ro 

■^ 

00 

<3" 

CN 

CN 

ro 

ro 

ro 

ro 

ro 

•^r 

rH 

0 

CN 

en 

-: 

J 

rtj 

00 

CO 

r-« 

CN 

rH 

rH 

rH 

rH 

H 

rH 

r- 

CN 

^0 

m 

en 

CQ 

IT) 

ifi 

rH 

en 

en 

rH 

in 

rH 

< 

en 

"^ 

00 

E-" 

CN 

rH 

«=r 

^P 

ro 

CO 

CN 

m 

CN 

in 

ro 

^3 

ro 

0 

CN 

CM 

CN 

rH 

CN 

CO 

rH 

en 

•^r 

ro 

ro 

CN 

in 

CM 

en 

VO 

0 

CQ 

r» 

r» 

<& 

r- 

r~ 

in 

^r 

p- 

in 

CN 

^0 

<-^ 

V£> 

^O 

en 

en 

en 

en 

cn 

en 

00 

"=r 

\Q 

0 

<o 

VD 

V£> 

^0 

UD 

vr> 

vc 

^0 

^0 

^0 

^0 

rH 

rH 

rH 

CN 

rH 

en 

en 

en 

en 

en 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

r-{ 

rH 

rH 

rH 

Q. 

<£> 

•rr 

^r 

ro 

r-\ 

en 

m 

^ 

CN 

^ 

r** 

^D 

m 

ro 

n-{ 

CN 

0 

m 

O 

**r 

CO 

■*r 

en 

r* 

-S* 

en 

CO 

0 

m 

CN 

< 

en 

en 

■^ 

r- 

00 

00 

ro 

rH 

00 

rH 

en 

l£ 

m 

CO 

r^ 

rH 

en 

en 

en 

r~- 

ro 

CN 

CN 

■^ 

rH 

^0 

ro 

^ 

^ 

r* 

r» 

vo 

^0 

<sO 

CN 

rH 

<-\ 

rH 

<-\ 

rH 

en 

^r 

^ 

00 

^r 

CN 

r-\ 

rH 

<-i 

r-{ 

rH 

rH 

rH 

rH 

<-\ 

n-\ 

•P 

H  C 

H 

(H 

H 

> 

> 

H 

M 

H 

> 

> 

H 

H 

H 

> 

> 

+ 

■r 

t 

J  -rH 

H  0 
:  Cm 

H 

H 
H 

H 

M 

H 
H 

H 

H 

H 
H 

H 

h 

H 

1+ 
( 

3 

5  C 

3 

t  PCS 

ro 

■^r 

in 

56 


ON 

CN 

■^ 

•** 

CN 

rH 

H 

in 

CN 

rH 

00 

T 

m 

•<* 

m 

Q) 

00 

00 

00 

00 

00 

CN 

CN 

CN 

CN 

CM 

CN 

CTi 

rH 

on 

en 

3 

en 

en 

CTt 

CTi 

en 

CN 

CN 

CN 

CN 

CN 

O 

en 

o 

en 

en 

H 

to 
> 

cq 

r» 

r- 

r** 

r* 

r-« 

H 

H 

rH 

H 

H 

r*» 

VO 

r» 

VO 

vo 

CTl 

cn 

en 

en 

en 

in 

in 

m 

m 

in 

en 

en 

en 

en 

en 

> 

•H 

ID 

o 

o 

m 

^ 

r^ 

H 

00 

vo 

en 

rr 

in 

^ 

CN 

r-« 

+J 

CT« 

CN 

o 

iH 

00 

•^ 

in 

"vT 

in 

r^ 

CN 

m 

H 

H 

00 

O 

< 

r>- 

O 

VO 

^ 

oo 

00 

CN 

CN 

CN 

CO 

H 

en 

o 

O 

00 

0) 

CN 

00 

O 

00 

<3< 

«^ 

H 

H 

rH 

00 

r^ 

vo 

r- 

r» 

00 

•n 

A 

H 

en 

H 

en 

o 

00 

m 

in 

in 

en 

en 

o 

o 

m 

CN 

O 

H 

H 

H 

H 

m 

r- 

VO 

en 

en 

rH 

VO 

CN 

oo 

in 

H 

r* 

r^ 

H 

r*» 

CQ 

TT 

CN 

O 

o 

oo 

O 

o 

<3« 

00 

o 

<3< 

VO 

oo 

m 

ro 

en 

en 

en 

CT> 

en 

00 

00 

P> 

r*. 

00 

^ 

in 

vo 

in 

rH 

LO 

in 

in 

in 

in 

vo 

VO 

VO 

vo 

VO 

00 

00 

CO 

00 

en 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

CN 

vo 

Q. 

ro 

p* 

CN 

ro 

CTt 

CN 

oo 

CN 

r>- 

^ 

CO 

00 

00 

r^ 

H 

r>« 

m 

vo 

en 

CN 

O 

r^ 

OO 

CTt 

^ 

in 

vo 

ro 

en 

r- 

< 

in 

p* 

vo 

CN 

H 

in 

m 

VO 

in 

oo 

in 

m 

vo 

in 

^r 

oo 

m 

n 

in 

^ 

00 

VO 

vo 

vo 

oo 

00 

00 

CO 

00 

CN 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

■^r 

^ 

00 

o 

oo 

T 

00 

00 

in 

o 

CO 

CO 

CN 

00 

00 

in 

r^ 

r- 

00 

00 

r» 

o 

o 

H 

rH 

o 

CN 

oo 

00 

rH 

CO 

m 

m 

m 

in 

m 

m 

•^ 

^« 

•"^ 

"5T 

<* 

r- 

vo 

m 

m 

VO 

o 

(N 

H 

in 

vo 

o 

vo 

in 

a\ 

•^ 

oo 

H 

en 

in 

VO 

CN 

^ 

< 

m 

rH 

r- 

m 

H 

in 

CN 

CM 

oo 

VO 

H 

CO 

oo 

H 

en 

in 

VO 

■^ 

r- 

ffl 

00 

«* 

"^ 

■^r 

CN 

r^ 

CN 

r*« 

en 

H 

vo 

^« 

H 

in 

CN 

en 

VO 

m 

m 

H 

oo 

m 

•^ 

O 

H 

H 

CO 

H 

m 

"5J1 

en 

VO 

OO 

o 

00 

VO 

ro 

H 

CN 

CN 

O 

"3* 

■** 

CO 

CN 

in 

H 

H 

o 

pq 

Ifl 

oo 

H 

H 

in 

■^r 

in 

■^ 

H 

■^ 

vo 

o 

en 

CTi 

VO 

o 

o 

o 

O 

o 

in 

m 

•^ 

m 

m 

00 

in 

m 

■q* 

r^ 

o 

O 

o 

O 

o 

00 

00 

00 

oo 

00 

ro 

oo 

ro 

00 

H 

rH 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

en 

iH 

en 

Q. 

r^ 

r- 

•*r 

C"\ 

CTl 

^ 

m 

CN 

o 

<sj< 

O 

CTl 

vo 

en 

in 

H 

in 

H 

00 

VO 

in 

r- 

CTl 

m 

CN 

r* 

oo 

CN 

vo 

vo 

i4 

CN 

o 

ro 

VO 

in 

oo 

ro 

CO 

en 

rH 

VO 

H 

en 

in 

vo 

CN 

00 

ro 

CN 

00 

VO 

CN 

CN 

H 

•^ 

CN 

m 

in 

in 

VO 

r- 

CT. 

r» 

en 

r» 

CO 

00 
H 

oo 

H 

oo 

rH 

00 

oo 
H 

00 
H 

00 

rH 

00 
rH 

CO 

in 

r™ 

^ 

R 

5  -P 

■r 

H    C 

H 

H 

H 

> 

> 

H 

H 

H 

> 

> 

H 

H 

H 

> 

> 

4- 

J  -H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

•r 

H    0 

H 

H 

H 

c 

:  cu 

i- 

i 

4- 

-i  en 

c 

)  c 

3 

VO 

r» 

00 

« 

■  ti 

57 


<D 
3 
H 

> 

> 

•H 
+J 
U 

<u 

•r-i 

o 


CQ 


cn 

CTi 

in 


vo 


VO 
■<* 
VO 


LD 

IT) 

LD 


in  ld 


ir»  in 

rH   rH 


VO 

m 
in 

m 


oo 

VO 

in 
in 


in 

m 


cn  oo  CO 

^  vo  vo 

r»  in  o 

m  in  oo 


in  in  Ln  in 


VO 


in 
cn 


rH 


m 
m 


CN 

rH 

m 
cn 


m 


m 


vo  CN 

^r   rH 

vo  in 


as 


o 

vo 


vo  in 


CN 

in 
cn 


m 

CO 

in 


cn     cn     n  m  m 


CQ 


Q. 


O 
CN 
CN 


CN 


oo  r» 

m  cn 

CN  rH 

CN  CN 


cn 

<3> 


oo  cn  m  oo 
p»  r^  r^  vo 


CN 

o  o 


CN   CN   CN   CN 
O   O   O   O 


00 
CO 


CN 


n 
o 

CN 
CN 


00 

CN 
CN 


CN 
O 
CN 
CN 


in 

CN 

cn 


<Ti 


vo 

rH 
CN 
CN 


^«  r-  i— l  cn  vo 

in  vo  c^o  oo  in 

cn  i-h  cn  rH  cn 

CN  CN  CN  CN  CN 


cq 


CN 

Q. 


in  cn  r* 
cn  co  cn 


o 


cn 
co 


CO   00 

<=r  o 


cn  cn 
vo  in 


oo 

vo 


CN   CN 

00   CO 


m 


l/l   (N   ^ 

o  cn  co 

CN   rH   rH 


cn  t~*  r-  vo 
vo  cn  cn  r-» 


o  o 
o  o 


cn 
cn 


o 


cn      rH 

a\     oo 


vo  r- 
m  oo 


o 
o 

CN 


oo 


o 


CN 

o 


CN 

cn 


o 


CQ 


^  tH  co  in 

CN  1"^  rH  co 

co  r~  cn  vo 

cn  cn  o  cn 


vo 
cn 


vo  cn  cn  <-\  cn 

o  co  ^  cn  ^ 

in  cn  rr  cn  o 

cn  cn  cn  cn  o 


00 


vo  vo  vo  vo 


r* 

vo 

r^ 

00 

r*» 

■^r 

CN 

00 

o 

<* 

in 

00 

in 

r^ 

CN 

cn 

cn 

o 

cn 

CO 

CN  CN  ^  Cn  rH 

r^  in  cn  oo  cn 

o  cn  "^  in  vo 

CN  CT>  rH  cn  O 


00 


cn 


vo 


vo 


id  +j 

•h  a 

-U  -H 
•H  0 

C  &4 


m     y-\     \-\     >     > 

H   H   H 


H   H   >   > 
\-4       M       l-i 


M-l  05 

o  a 
3 

*  cd 


cn 


58 


TABLE  IV 
Actual  values  X   =  0.1,  a  =  1.0,  6  =  0.05 


#  of 
Runs 

Initial- 
Point 

A 

a 

3 

A 

B 

A 

B 

A 

B 

1 

I 

II 

III 

IV 

V 

0.1202 
0.1199 
0.1168 
0.1158 
0.1203 

0.1203 
0.1203 
0.1203 
0.1204 
0.1201 

343.512 

129.591 

8.898 

6.577 

4.182 

676.800 
968.150 
1506.926 
874.957 
955.629 

0.0400 
0.0401 
0.0409 
0.0412 
0.0400 

0.0400 
0.0401 
0.0400 
0.0385 

0.040.1 

2 

I 

II 

III 

IV 

V 

0.0984 
0.0865 
0.0977 
0.0798 
0.0500 

0.0863 
0.0863 
0.0  86  3 
0.0863 
0.0863 

7.534 
1.306 
6.158 
0.810 
0.154 

1.297 
1.296 
1.297 
1.292 
1.296 

0.0478 
0.0509 
0.0480 
0.0527 
0.0583 

0.0510 
0.0510 
0.0510 
0.0511 
0.0511 

3 

I 

II 

III 

IV 
V 

0.1390 
0.1430 
0.1440 
0.1434 
0.0438 

0.1434 
0.1434 
0.1436 
0.1433 
0.1434 

8.123 
138.531 
204.078 
414.903 

0.055 

392.550 
1520.067 
1465.586 
1177.650 
2056.280 

0.0343 
0.0332 
0.0328 
0.0331 
0.0515 

0.0331 
0.0332 
0.0330 
0.0  332 
0.0331 

4 

I 

II 

III 

IV 
V 

0.0878 
0.0890 
0.0891 
0.0870 
0.0899 

0.0843 
0.0873 
0.0856 
0.0829 
0.0860 

0.118 
0.12  3 
0.12  3 
0.115 
0.125 

0.10  7 
0.116 
0.111 
0.10  3 
0.112 

0.0394 
0.0392 
0.0393 
0.0  39  4 
0.0390 

0.0397 
0.0393 
0.0395 
0.0398 
0.0394 

5 

I 

II 
III 

IV 
V 

0.0507 
0.0695 
0.0692 
0.0540 
0.0676 

0.0508 
0.0509 
0.0508 
0.0508 
0.0507 

4.944 
0.847 
3.882 
0.826 
0.672 

0.509 
0.513 
0.508 
0.509 
0.506 

0.0530 
0.0988 
0.0749 
0.0717 
0.10  34 

0.0775 
0.0775 
0.0775 
0.9775 
0.0776 
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TABLE 

IV   Cont. 

#    of 
Runs 

Initial 
Point 

A 

ex 

3 

A 

B 

A 

B 

A 

B 

I 

0.1384 

0.0994 

9.0  72 

0.570 

0.0494 

0.0591 

II 

0.1019 

0.0996 

0.628 

0.576 

0.0585 

0.0590 

6 

III 

0.1363 

0.0998 

6.105 

0.579 

0.0499 

0.0589 

IV 

0.10  79 

0.0998 

0.810 

0.582 

0.0571 

0.0590 

V 

0.1272 

0.0994 

2.439 

0.571 

0.0525 

0.0590 

I 

0.115  8 

0.0738 

6.780 

0.301 

0.0405 

0.0502 

II 

0.0755 

0.0738 

0.324 

0.301 

0.0496 

0.0502 

7 

III 

0.0752 

0.0743 

0.319 

0.309 

0.0499 

0.0501 

IV 

0.0  75  7 

0.0739 

0.328 

0.30  8 

0.0499 

0.0501 

V 

0.1188 

0.0738 

28.839 

0.301 

0.0397 

0.0502 

I 

0.0753 

0.0748 

7.325 

5.783 

0.0645 

0.0631 

II 

0.0739 

0.0740 

4.647 

4.763 

0.0634 

0.0634 

8 

III 

0.0735 

0.0735 

4.219 

4.221 

0.0635 

0.0635 

IV 

0.0737 

0.0741 

4.361 

4.831 

0.0634 

0.0633 

V 

0.0170 

0.0010 

0.033 

0.0008 

0.0724 

0.0634 

I 

0.1257 

0.1252 

11.646 

9.302 

0.0275 

0.0276 

II 

0.1252 

0.125  3 

10.181 

9.374 

0.0275 

0.0275 

9 

III 

0.1240 

0.1246 

6.983 

8.598 

0.0278 

0.0278 

IV 

0.1254 

0.1254 

11.010 

10.296 

0.0276 

0.0275 

V 

0.1017 

0.1255 

0.754 

10.407 

0.0330 

0.0275 

I 

0.1387 

0.1438 

5.688 

25.1.70 

0.0326 

0.0313 

II 

0.1441 

0.1428 

28.830 

29.290 

0.0312 

0.0310 

10 

III 

0.1400 

0.1440 

6.290 

27.747 

0.0326 

0.0313 

IV 

0.1437 

0.1441 

21.370 

26.590 

0.0313 

0.0312 

V 

0.1415 

0.1427 

10.0  30 

14.610 

0.0310 

0.0316 
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C.   TEST  PROCEDURE  FOR  PARAMETERS 

In  order  to  determine  whether  the  model  reasonably  fits 
the  data  a  simple  test  procedure  is  applied.   It  is  to 
exhibit  the  simple  comparison  of  the  actual  points   z(P) 
and  the  estimated  values   z(P),  using  the  estimated 
parameters . 

Basically   z (P)   can  be  obtained  as  follows: 

x(p)  =  -  J  An(l-p)  (4.15) 

Then  if  it  is  substituted  in  expression  (2.37),  z (P)   will 
be 

z(P)  = ox_JP) (4.16) 

(1  +  ax(P)  )  (1  +  J3x(P)  ) 

In  equation  (4.16) ,  all  variables  are  known  so  estimated 
quantiles  can  be  easily  obtained. 

In  Fig.  9  actual  values  and  estimated  values  are  plotted, 
Also  Fig.  10  displays  residuals,  which  are  the  difference 
between  actual  and  estimated  values. 
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V.   NUMERICAL  APPLICATIONS  TO  TWO  SETS 
OF  REAL  DATA 

In  this  section  the  two  failure  data  sets  studied  come 
from  these  sources:   Oral  irrigator  [7]  and  human  life  [8]. 
They  are  used  for  numerical  examples  in  that  these  data  are 
fitted  using  the  model  (2.6). 

A.   ORAL  IRRIGATOR 

The  data  used  was  obtained  from  the  Commun.  Statist. - 
Theor.  meth. ,  Colvert  and  Bordman  [7].   The  data  was  col- 
lected such  that  100  oral  irrigators  were  placed  on  the  test 
was  terminated  at  700  time  units.   During  the  life  test,  9  8 
oral  irrigators  failed;  another  2  oral  irrigators  survived. 
The  ordered  observed  times  to  failure  of  the  oral  irrigators 
is  tabulated  in  Table  V. 

Using  the  data  of  Table  V,  the  parameters   a,  8f  and  A 
are  estimated  by  means  of  the  3-percentile  approach.   To 
demonstrate  the  differences  between  estimated  values, 
various  different  combinations  of  p  values  were  used.  The 
results  are  shown  in  Table  VI. 

Examination  of  Table  VI  indicates  that  the  3-percentile 
approach  may  produce  estimates  having  great  differences  for 
different  percentile  values,  except  here  in  the  case  of  the 
first  two  combinations.  Also  the  first  two  estimations  seem 
to  have  accepatable  limiting  age.  Moreoever,  the  last  two 
combinations  do  not  fit  well,  as  judged  from  the  residuals. 
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TABLE  V 
TIME  TO  FAILURE  OF  ORAL  IRRIGATORS 


1.75 

7.02 

7.58 

9.76 

15.02 

15.5  7 

17.39 

19.55 

22.47 

23.24 

23.96 

25.05 

32.44 

36.87 

42.76 

43.14 

46.95 

56.33 

58.99 

59.08 

60.37 

61.01 

•    77.86 

86.45 

88.50 

103.06 

104.34 

105.85 

117.46 

120.11 

122.28 

122.61 

129.31 

130.42 

137.57 

142.27 

142.98 

148.29 

150.79 

151.21 

155.62 

157.93 

160.72 

169.79 

186.26 

197.60 

224.83 

233.64 

242.07 

256.86 

260.77 

261.68 

277.99 

283.95 

288.94 

295.48 

314.76 

316.06 

332.07 

339.46 

362.61 

369.47 

370.74 

491.06 

403.39 

414.78 

426.71 

459.62 

455.84 

457.94 

466.61 

468.64 

469.09 

476.42 

481.41 

481.82 

488.15 

490.06 

493.67 

494.38 

503.72 

508.93 

509.01 

418.32 

532.29 

534.62 

545.23 

547.41 

558.41 

571.10 

585.52 

539.11 

592.93 

607.15 

623.15 

647.91 

TABLE  VI 
ESTIMATED  VALUES  OF  a,    &,    \      BY  USING  3-PERCENTILE 


Pi 

P2 

P3 

X 

a 

e 

1/6 

.1 

.25 

.1 
.25 

.5 
.5 

.5 
.5 

.9 
.9 

.75 
.75 

0.00153 
0.00158 
0.00239 
0.00273 

0.00634 
0.00710 
0.01835 
0.04514 

0.001028 
0.001018 
0.000228 
0.000077 

972.76 

982.31 

4385.96 

12987.01 
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Next  the  nonlinear  optimization  approach  was  applied  to 
estimate  parameters   a,  3   and  A   associated  with   p  ,   p„ 
and   p^,  both  using  the  derivative  and  without  derivative  by 
taking  into  account  9  8  ordered  values.   The  results  are 
tabulated  in  Tables   VII  and  VIII  for   p,,   p2,   p    and  X 
a,      3  respectively.   Some  initial  points  in  Table  II  are 
used.   Table  VII  indicates  that  the  estimates  obtained  from 
the  nonlinear  estimation  method  are  more  consistent  than 
those  from  the  3-percentile  approach.   Also,  the  boundary 
points ,  1/3/  seem  reasonable.   Finally,  comparison  of  the 
actual  values  of   z (p)   and  estimated  values  of   z(p)  indi- 
cate that  the  fitting  is  reasonable.    The  residuals  of 
fitted  model  and  the  comparison  of  the  actual  and  estimated 
values  are  plotted  in  Fig.  11  and  Fig.  12. 

B.   HUMAN  LIFE  (MORTALITY)  DATA 

The  data  used  in  this  example  was  obtained  from  the 
1969-71  life  table  [8]  for  white  females  in  the  United  States 
The  table  has  been  prepared  from  a  history  of  100,000  persons 
the  number  of  surviving,  and  the  number  dying  has  been  given 
for  each  age  interval.   Look  at  Table  IX. 

1.   Estimation  by  Using  3-Percentile  Approach 

The  life  data  is  first  fitted  to  the  model  (2.37) 
by  using  the  3-percentile  approach.   Certain  pth  percentile 
values  are  selected  and  used  in  the  estimation  process. 
These  results  are  tabulated  in  Table  X.   Investigation  of 
the  last  eight  combinations  of  percentiles  in  these  tables 
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indicates  that  the  limiting  age,  1/S ,  is  approximately  110, 
which  is  reasonable.   However,  the  first  combination  of 
percentiles  gives  a  good  fit  in  terms  of  the  sum  of  the 
squared  errors.   But  for  these  limiting  age  is  reduced  to 
about  95  years.   This  apparently  means  that  the  fit  of  the 
model  does  not  represent  the  age  above  95.   Thus  it  can 
be  said  that  there  is  a  trade-off  between   1/3   and  the  sum 
of  squared  errors.   The  present  model  simply  does  not  seem 
to  fit  the  mortality  data  very  well. 

2 .   Estimation  By  Using  Nonlinear  Least  Squares,  and 
Application  of  a  Nonlinear  Programming  Algorithm 
In  the  previous  section  we  discussed  the  fact  that 
one  of  the  problems  encountered  in  the  3-percentile  estima- 
tion procedure  was  the  high  uncertainty  of  fitting,  as 
measure  by  the  sum  of  squared  errors.   In  order  to  reduce 
this  variability,  a  nonlinear  estimation  process  is  applied. 
Previously  it  was  noted  that  the  GRG  package  can  be  used 
with  analytically  computed  derivatives  and  without 
derivatives.   If  it  is  used  with  analytically  computed 
derivatives,  it  is  necessary  to  provide  another  subprogram 
by  the  user.   Otherwise  the  GRG  package  will  provide  the 
derivative,  computed  numerically  and  automatically. 

The  derivative  of  the  objective  function  (4.13)  is 
the  normal  equations  such  that 
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TABLE    IX 


M 

ORTALITY    DAT 

A  FOR  WHITE    FEMALES 

,    1969-71 

Time 

#        Of 

#    Of 

Time 

#   of 

#    of] 

Interval 

Survivors 

Dying 

Interval 

Survivors 

Dying 

0-    1 

100,000 

15  32 

27-28 

97,165 

70 

1-   2 

99,468 

100 

28-29 

97,095 

73 

2-    3 

98,368 

65 

29-30 

97,022 

77 

3-   4 

98,303 

54 

30-31 

96,945 

81 

4-   5 

98,249 

46 

31-32 

96,864 

87 

5-   6 

98,203 

39 

32-33 

96,777 

93 

6-    7 

98,164 

35 

33-34 

96,684 

101 

7-    8 

98,129 

32 

34-35 

96,583 

109 

8-    9 

98,097 

29 

35-36 

96,474 

118 

9-10 

98,068 

26 

36-37 

96,356 

128 

10-11 

98,042 

24 

37-38 

96,228 

141 

11-12 

98,018 

24 

38-39 

96,087 

153 

12-13 

97,994 

25 

39-40 

95,932 

170 

13-14 

97,969 

30 

40-41 

95,762 

185 

14-15 

97,939 

37 

41-42 

95,577 

201 

15-16 

97,902 

45 

42-43 

95,376 

220 

16-17 

97,857 

54 

43-44 

95,156 

242 

17-18 

97,803 

60 

44-45 

94,914 

263 

18-19 

97,743 

62 

45-46 

94,649 

291 

19-20 

97,681 

63 

46-47 

94,358 

317 

20-21 

97,618 

63 

47-48 

94,041 

344 

21-22 

97,555 

63 

48-49 

93,697 

372 

22-23 

97,492 

63 

49-50 

93,325 

401 

23-24 

97,429 

64 

50-51 

92,924 

433 

24-25 

97,365 

66 

51-52 

92,491 

469 

25-26 

97,299 

66 

52-53 

92,022 

506 

26-27 

97,233 

68 

53-54 

91,516 

546 
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TABLE    IX   Cont. 


Time 

#    of 

#    of 

Time 

#   of 

#    of 

Interval 

Survivors 

Dying 

Interval 

Survivors 

Dying 

54-55 

90,970 

587 

82-    83 

41,215 

3,586 

55-56 

90,383 

632 

83-    84 

32,629 

3,589 

56-57 

89,751 

680 

84-    85 

34,040 

3,550 

5  7-5  8 

89,071 

730 

85-    86 

30,490 

3,495 

58-59 

88,341 

781 

86-    87 

26,995 

3,425 

59-60 

87,560 

834 

87-    88 

23,570 

3,286 

60-61 

86,726 

911 

88-    89 

20,284 

3,0  72 

61-62 

85,835 

953 

89-    90 

17,212 

2,806 

62-63 

84,882 

1,022 

90-    91 

14,406 

2,531 

63-64 

83,860 

1,098 

91-   92 

11,875 

2,262 

64-65 

82,762 

1,183 

92-   93 

9,613 

1,982 

65-66 

81,579 

1,275 

93-    94 

7,63.1 

1,694 

66-67 

80,304 

1,375 

94-    95 

5,937 

1,411 

67-68 

78,929 

1,486 

95-    96 

4,526 

1,145 

68-69 

77,443 

1,607 

96-    97 

3,381 

905 

69-70 

75,836 

1,735 

97-    98 

2,476 

696 

70-71 

74,101 

1,862 

98-    99 

1,780 

524 

71-72 

72,239 

1,993 

99-100 

1,256 

384 

72-73 

70,246 

2,141 

100-101 

872 

277 

73-74 

68,105 

2,312 

101-102 

595 

195 

74-75 

65,793 

2,503 

102-103 

400 

135 

75-76 

63,290 

2,693 

103-104 

265 

92 

76-77 

60,597 

2,872 

104-105 

173 

61 

77-78 

57,725 

3,0  39 

105-106 

112 

41 

78-79 

54,686 

3,187 

106-107 

71 

26 

79-80 

51,499 

3,317 

107-108 

45 

17 

80-81 

48,182 

3,434 

108-109 

28 

11 

81-82 

44,748 

3,533 

-i- 

109-110 

17 

6 
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TABLE  X 

ESTIMATED  VALUES  OF   a,   3,  A   FOR  MORTALITY  DATA  USING 

3 -PERCENTILE  APPROACH 


Percentile 

/K 

/% 

/\ 

Values 

X 

a 

3 

•1/3 

22      56 

92 

0.000674 

0.1190 

0.010550 

94.7 

20      56 

92 

0.000550 

0.0404 

0.010520 

95.0 

17      56 

92 

0.000508 

0.0302 

0.010550 

94.7 

16      56 

92 

0.000440 

0.0199 

0.010550 

94.7 

17      56 

109 

0.000862 

0.1676 

0.009069 

110.2 

18      56 

109 

0.000886 

0.2389 

0.009068 

110.2 

19      56 

109 

0.000908 

0.3708 

0.009067 

110.2 

17      56 

110 

0.000876 

0.1805 

0.008991 

111.2 

18      56 

110 

0.000899 

9.2613 

0.008989 

111.2 

19      56 

110 

0.000920 

0.4188 

0.008988 

111.2 

20      56 

110 

0.000941 

0.9070 

0.008988 

111.2 
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8S__  »    (  P2^(Pi) ) 

3Pl  i=l{(i+    P2£(Pi})  (1    +    P3e<Pi))  ( 


=    -2 


J  P1P2£2(P.)  ) 

x  j2(pi}  -  (i  +  p2e(p.))(i  +  p3s(Pi))j       <4-17> 


N      (Pjg     (Pj)  (1    +    P2£(Pi))     -£(Pi)     PlP2£  ~(PJ| 
8p2  i=l  I  (1    +   p2£(Pi))2     (1    +   p3e(Pi))  J 

j  PXP2£     (P.)  ) 

X    (Z(Pi3    "    (1    +p2£(P.))(l1+P3£(P,)))  (4«18) 


3^=    +    2      j      j  PjP^V | 


3 


=  1    (     (1     +     P2£(Pi))2       (1     +     p3£(Pi))     j 


j  PlP     £2(P     )  ) 

X   j    Z(Pi)     "     (1    +    p2£(P.))(l    +    P3£(Pi))j  (4'19) 

N   is  110  in  the  case  of  human  mortality  data. 

The  results  of  estimation  were  shown  in  Tables  XI 
and  XII  for   p.   and  a,   g,  A   respectively.   Also  Fig.  13 
and  Fig.  14  demonstrate  the  actual  and  estimated  values  and 
residuals  of  fitting. 

Examination  of  Table  XI  and  Table  XII  indicates 
that  the  objective  function  values  (the  sum  of  the  squared 
errors)  of  the  fitting  are  smaller  than  the  objection  func- 
tion values  of  3-percentile  approach.   Also  the  estimated 
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parameter  values   3   and  X      associated  with   p    and  p 
are  very  close  to  each  other.   But  estimated  values  of  a 
show  considerable  difference  for  the  various  initial  values. 
However,  there  are  two  important  facts  to  notice.   First, 
when  a   increases  significantly,  the  objective  function 
values  remain  almost  constant.   That  is,  it  is  not  effected 
significantly  on  the  objective  function  values  in  this  case. 
Second,  this  analytical  model  (2.6)  does  not  represent 
deaths  beyond  the  age  95..  However,  modifications  can  be 
made  in  the  model  (2.6)  for  this  kind  of  difficulty.   Some 
ideas  will  be  discussed  in  a  later  section. 

3 .   Model  Modifications 

In  the  previous  section  it  was  noted  that  the  model 
(2.6)  has  relatively  great  errors  and  does  not  accurately 
describe  probability  of  death  at  age  greater  than  95  in  the 
human  life  example.   This  means  that  the  hazard  function 
h  (z)   increases  too  rapidly  in  the  wearout  period.   If  it 
is  possible  to  slow  the  rate  of  increase  of  hazard  perhaps  a 
better  result  can  be  obtained. 

In  Section  II,  it  was  stated  that  the  function 
R(x)   describes  the  wearout  period  in  bath  tub  type  curves: 

R(x)  "  5-hi5  '      6  >  ° 

If  a  new  parameter,   y,  is  defined  which  is  between 
0  and  1  and   R1 (x)   is  now  defined  to  be  the  y    power  of 
R(x)  : 
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R'  (x)  =  ~,      g  >  o,   0  <  y  <  1  (4.20) 

(1  +  3x)y 


then   R1 (x)   will  provide  for  a  slowly  increasing,  rather 
than  rapidly  increasing,  in  wearout  period.   With  this 
revision  the  model  (2.6)  now  becomes 

z  =  G(x)  =  xL(x)-R'  (x)  (4.21) 

or 

2 

z  = (4.22) 

(1  +  ax)  (1  +  3x)y 

where   a  >  0,   3  >  0   and  0  <  y   <   1. 

However,  after  we  put  another  variable  in  the  model, 
it  will  be  too  difficult  to  handle  in  the  previous  estimation 
technique  for  obtaining  values  of  the  parameters   a,   3,  A 
and  y   because  of  nonlinearity  and  indef initeness  of  the 
expression  (4.22). 

For  this  reason,   y  will  be  assumed  constant  in  the 
previous  estimation  procedure.   Then  it  will  be  computa- 
tionally convenient.   Actually,  the  estimation  procedure 
does  not  require  any  change.   The  nonlinear  least  square 
estimation  approach  will  simply  be  used  for  different  values 
of  y*   Table  XIII  and  Fig.  15  demonstrate  the  parameters 
estimated  and  the  resulting  objective  function  values.   From 
Table  XIII,  it  can  be  easily  seen  that  the  new  result  has  a 
smaller  sum  of  squared  errors.   The  best  value  of   y   seems 
to  be  near  0.95. 
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TABLE  XIII 
ESTIMATION  VALUES  OF   a,   3,  X 
USING  CONSTANT  y 


#   of 

Objective 

Run 

Y 

A 

a 

3 

Value 

1 

1.000 

0.000834 

0.42950 

0.01029 

4,335.89 

2 

0.97 

0.000512 

0.04051 

0.01182 

3,119.85 

3 

0.95 

0.000363 

0.01928 

0.01314 

2,494.10 

4 

0.93 

0.000740 

1.06371 

0.01352 

2,999.11 

5 

0.90 

0.000685 

1.81283 

0.01554 

3,267.47 

6 

0.85 

0.000592 

5.74169 

0.02031 

5,099.37 

7 

0.7 

0.000319 

3.17467 

0.04580 

15,310.06 
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C.   CONCLUSIONS 

Simple  analytical  hazard  models  have  been  developed 
and  fitted  to  situations  (data)  that  exhibit  bath  tub 
shaped  hazard  functions.   That  is,  failure  rates  may  be 
high  at  early  ages  ("infant  mortality") ,  constant  at 
intermediate  ages,  and  high  again  for  later  ages  ("wearout") 
The  procedure  emphasizes  representations  of  the  inverse 
distribution  function;  simulation  is  thus  facilitated. 

The  failure  time  distributions  so  derived  should  be 
useful  in  analyzing  maintenance  and  replacement  policies. 

A  least  squares  technique  for  fitting  the  hazard  models 
to  data  are  suggested  and  applied. 
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